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ABSTRACT 

We have simulated the formation of a galaxy cluster in a ACDM universe using twelve differ¬ 
ent codes modeling only gravity and non-radiative hydrodynamics (ART, Arepo, Hydra 
and 9 incarnations of Gadget). This range of codes includes particle based, moving and 
fixed mesh codes as well as both Eulerian and Lagrangian fiuid schemes. The various Gad¬ 
get implementations span traditional and advanced smoothed-particle hydrodynamics (SPH) 
schemes. The goal of this comparison is to assess the reliability of cosmological hydrodynam- 
ical simulations of clusters in the simplest astrophysically relevant case, that in which the gas 
is assumed to be non-radiative. We compare images of the cluster at ^ = 0, global properties 
such as mass, and radial profiles of various dynamical and thermodynamical quantities. The 
underlying gravitational framework can be aligned very accurately for all the codes allowing a 
detailed investigation of the differences that develop due to the various gas physics implemen¬ 
tations employed. As expected, the mesh-based codes ART and Arepo form extended en¬ 
tropy cores in the gas with rising central gas temperatures. Those codes employing traditional 
SPH schemes show falling entropy profiles all the way into the very centre with correspond¬ 
ingly rising density profiles and central temperature inversions. We show that methods with 
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1 INTRODUCTION 

Galaxy clusters are the largest virialized objects in the Universe 
and, as such, provide both an important testbed for our theories 
of cosmological structure formation and a challenging laboratory 
within which to study the fundamental physical processes that drive 
galaxy formation and evolution. The overdense regions that con¬ 
tain clusters in the present-day Universe were the first to collapse 
and virialize in the early Universe, and so our theories must pre¬ 
dict their assembly history over a large fraction of the lifetime of 
the Universe (see Allen, Evrard & Mantz 2011 and Kravtsov & 
Borgan! 2012 for a recent review). At the same time, galaxies in 
the cores of clusters have orbited within the often violently grow¬ 
ing cluster potential over many dynamical times, while the broader 
cluster galaxy population is continually replenished by the infall of 
galaxies from the field. 

Computer simulations are now well established as a powerful 
and indispensable tool in the interpretation of astronomical obser¬ 
vations (see for instance Borgani & Kravtsov 2011). Early A^-body 
simulations (White 1976; Eall 1978; Aarseth, Turner & Gott 1979), 
which included the gravitational effects of dark matter only, were 
vital in interpreting the results of galaxy redshift surveys and un¬ 
veiling the large scale structure of the Universe, and subsequently 
in resolving structure in the non-linear regime of dark matter halo 
formation. The focus of modern simulations has now shifted to 
modeling galaxy formation in a cosmological context (see Borgani 
& Kravtsov 2011 for a recent review), incorporating the key phys¬ 
ical processes that drive galaxy formation - such as the cooling of 
a collisional gaseous component (e.g. Pearce et al. 2000; Muan- 
wong et al. 2001; Dave, Katz & Weinberg 2002; Kay et al. 2004; 
Nagai, Vikhlinin & Kravtsov 2007; Wiersma, Schaye & Smith 
2009), the birth of stars from cool overdense gas (e.g. Springel 
& Hernquist 2003; Schaye & Dalla Vecchia 2008), the growth of 
black holes (Di Matteo, Springel & Hernquist, 2005), and the in¬ 
jection of energy into the inter-stellar medium by supernovae (e.g. 
Metzler & Evrard 1994; Borgani et al. 2004; Dave, Oppenheimer 
& Sivanandam 2008Dalla Vecchia & Schaye 2012) and powerful 
AGN outflows (e.g. Thacker, Scannapieco & Couchman 2006; Si- 
jacki et al. 2007; Puchwein, Sijacki & Springel 2008; Sijacki et al. 
2008; Booth & Schaye 2009). These processes span an enormous 
dynamic range, both spatial and temporal, from the sub-pc scales 
of black hole growth to the accretion of gas onto Mpc scales from 
the cosmic web. Galaxy clusters offer an ideal testbed for the study 
of these processes and their complex interplay, precisely because 
their enormous size encompasses the range of scales involved. Eor 
this reason, the study of galaxy formation and evolution in cluster 
environments occupies a fundamental position in observational and 
numerical astrophysics. 

This raises the important question of how reliably cosmolog¬ 
ical galaxy formation simulations recover the properties of galaxy 
clusters. In the now classic Santa Barbara Cluster Comparison 
Project, Erenk et al. (1999) (SB99 from now on) compared the 
bulk dark matter and gas properties of a galaxy cluster formed in 
a non-radiative cosmological hydrodynamical simulation run using 
twelve (then state-of-the-art) mesh- and particle-based (hereafter 
smoothed particle hydrodynamics, SPH) codes. They displayed vi¬ 
sual representations of the cluster at z = 0 and at z = 0.5 when 
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a major merger event was ongoing, and compared several observ¬ 
able quantities such as enclosed mass, gas temperature and X-ray 
emission. The large scatter in essentially all quantities between the 
twelve models is evident from the plots. The origin of these discrep¬ 
ancies was partly the poor timing agreement between the methods 
due to an ambiguity in the start redshift as well as large differences 
in the effective numerical resolution that arose because the sim¬ 
ulation challenged the computing resources available at the time. 
The key discrepancy has, however, stood the test of times: the di¬ 
vergence between mesh-based and SPH codes in the radial entropy 
profile of the gas, defined in SB99 as 

S{R) = log [Tgas(-R)/Pgas(-R)"^"] (1) 

where R is the spherical radius with respect to the cluster centre 
of mass; Tgas is the gas temperature; and pgas is the gas density. 
Pig. 18 of SB99 showed some tentative indication that the entropy 
profiles of the grid based codes (principally those of Bryan and 
Cen) displayed a constant entropy core whereas the entropy profiles 
of the sparticle based SPH codes continued to fall all the way into 
the centre. 

These results have been confirmed subsequently by several 
studies (Voit, Kay & Bryan, 2005; O’Shea et al., 2005; Dolag et al., 
2005; Wadsley, Veeravalli & Couchman, 2008; Mitchell et al., 
2009). Eor example, Wadsley, Veeravalli & Couchman (2008) and 
Mitchell et al. (2009) suggested that the discrepancy owes also 
to the artificial surface tension and the associated lack of multi¬ 
phase fluid mixing in classic SPH, while similar conclusions have 
been reached by Sijacki et al. (2012) when comparing the moving 
mesh code Arepo of Springel (2010) with classic SPH, using P- 
GADGET3 with the entropy conserving SPH version of Springel 
& Hernquist (2002). Interestingly, in their recent study. Power, 
Read & Hobbs (2014) have shown that SPHS (Read & Hayfield, 
2012a), an extension of SPH that improves the treatment of mix¬ 
ing by means of entropy dissipation, produces constant entropy 
cores in non-radiative galaxy cluster simulations that are consis¬ 
tent with the results of the adaptive mesh refinement (AMR) code 
RAMSES (Teyssier, 2002). Both Wadsley, Veeravalli & Couch¬ 
man (2008) and Maier et al. (2009) report entropy cores in runs that 
incorporate sub-grid models for turbulence. These results suggest 
that modern SPH codes can overcome the problems first reported 
in Erenk et al. (1999) and subsequently in Agertz et al. (2007). It 
is worth noting that it is not obvious that mesh-based codes neces¬ 
sarily recover the correct form for the entropy profile. Eor example, 
Springel (2010) reports significant variation in the entropy profile 
for the same AMR code (ENZO) that is particularly sensitive to 
choice of refinement criteria. More generally, differences are ap¬ 
parent when comparing AMR results to that of the moving mesh 
code Arepo Springel (2010) report an entropy core that is signif¬ 
icantly lower than that found in AMR codes (e.g. compare Eigure 
45 of Springel 2010 with Eigure 18 of Erenk et al. 1999 or Eigure 
5 of Voit, Kay & Bryan 2005). 

In this work - emerging out of the ‘nIETy cosmology’ work¬ 
shop^ - we revisit the idea of the Santa Barbara Cluster Compar¬ 
ison Project fifteen years later. We take four modern cosmolog¬ 
ical simulation codes (with one of them taken in nine different 
flavors, for a total of twelve different codes) and study the for¬ 
mation and evolution of a large galaxy cluster (with final mass 
^200 — 1.1 X lO^^M©). Eirst we perform blind dark matter only 
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Table 1. List of all the simulation codes participating in the nIFTy cluster 
comparison project. 


Code name 

Reference 

CART 

Rudd, Zentner & Kravtsov (2008) 

Arepo 

Springel (2010) 

Hydra 

Couchman, Thomas & Pearce (1995) 

Gadget: 

Springel (2005) 

G2-Anarchy 

Dalla Vecchia et al. in prep. < 

G3-X 

Beck et al. (2015) 

G3-SPHS 

Read & Hayfield (2012a) 

G3-Magneticum 

Hirschmann et al. (2014) 

G3-PESPH 

Huang et al. in prep. 

G3-MUSIC 

Sembolini et al. (2013) 

G3-OWLS 

Schaye et al. (2010) 

G2-X 

Pike et al. (2014) 


simulations with the favored parameter sets of each group. The re¬ 
sults from these show the typical scatter that is expected for cur¬ 
rently published works in this area. We then use a common param¬ 
eter set to align our gravitational framework before continuing to 
study non-radiative gas simulations. This allows us to focus solely 
on the differences between the models that arise due to the different 
hydrodynamical implementations. This also permits us to cleanly 
study the formation (or not) of a gas entropy core. 

The rest of this paper is organized as follows: in Section 2 we 
briefly describe the twelve methods used in this study and supply 
tables listing parameter choices. In Section 3 we describe how our 
initial conditions were generated and chosen. The main results are 
presented in Section 4, which discusses the dark matter only simu¬ 
lations, and in Section 5, which presents the results from the non- 
radiative simulations. We discuss convergence and scatter among 
the different codes in Section 6. We summarize out results in Sec¬ 
tion 7. We present additional supporting material the Appendix. 


2 THE CODES 

The numerical codes used for this project can be divided into two 
main groups: mesh-based and SPH codes. The mesh based codes 
used in this work are ART (Kravtsov, Klypin & Khokhlov, 1997) 
and Arepo (Springel, 2010): the first one uses Eulerian hydrody¬ 
namics as the second one uses a moving unstructured mesh and 
can be considered almost Lagrangian. SPH codes use Lagrangian 
hydrodynamics: we used Hydra (Couchman, Thomas & Pearce, 
1995) and nine different versions of Gadget (Springel, 2005). 
Among SPH codes we can distinguish classic and modern SPH, 
defining the latter as codes that adopted an improved treatment 
of discontinuities. The codes employ different techniques to solve 
the evolution equations for a two-component fluid of dark matter 
and non-radiative gas coupled by gravity. To calculate gravitational 
forces, ART uses Adaptive Mesh Refinement (AMR), Arepo and 
Gadget are based on TreePM (Tree-i-Particle-Mesh) methods and 
Hydra uses AP^M (Adaptive Particle-Particle, Particle-Mesh). 
Gas particles are treated in the following ways: by means of SPH in 
Gadget and Hydra , using a Voronoi mesh in Arepo, and using 
Eulerian AMR in ART. 

The following short summaries detail the specifics of each 
simulation code contributing to this comparison (the reference au¬ 
thor for each code is the person who run the simulation). We focus 


on key differences and novel aspects. Generalized descriptions of 
the individual codes can be found in their respective methods pa¬ 
pers. Table 1 provides the standard reference for each code; Table 
2 summarizes the key numerical parameters used for each run. Sec. 
2.3 describes the main improvements made by modern SPH codes. 


2.1 Mesh-based codes 


CART (Nagai, Nelson & Lau) 

ART (Adaptive Refinement Tree, ART) is an N-body plus 
hydrodynamics code with adaptive mesh refinement. The code 
solves the inviscid fluid dynamical equations using a second or¬ 
der accurate Godunov method with piecewise-linear reconstructed 
boundary states and the exact Riemann solver of Colella & Glaz 
(1985). The current version of the code used for this work is CART 
(Chicago-ART), which it is parallelized for distributed machines 
using MPI and features a flexible time-stepping hierarchy. 

Arepo (Puchwein) 

Arepo uses a Godunov scheme on an unstructured moving 
Voronoi mesh. The mesh cells move (roughly) with the fluid. The 
main differences to Eulerian AMR codes consist in that AREPO 
is almost Lagrangian and it is Galilean invariant by construction; 
furthermore, AREPO has automatic refinement for hydrodynamics 
and gravity and uses a Tree-PM gravity solver. The main difference 
to SPH codes is that the hydrodynamic equations are solved with 
a finite-volume Godunov scheme. In this work, we always use the 
total energy as a conserved quantity in the Godunov scheme rather 
than the entropy-energy formalism described in Springel (2010). 


2.2 SPH codes 


Gadget2-Anarchy (Dalla Vecchia) 

Gadget-Anarchy (G2-Anarchy) is an implementation of Gad- 
GET2 employing the pressure-entropy SPH formulation derived by 
Hopkins (2013). The artificial viscosity switch has been been im¬ 
plemented following Cullen & Dehnen (2010), whose algorithm 
allows precise detection of shocks and avoid excessive viscosity in 
pure shear flows. G2-Anarchy uses a purely numerical switch for 
entropy diffusion similar to the one of Price (2008), but without 
requiring any diffusion limiter. The kernel adopted is the func¬ 
tion of Wendland (1995) with 100 neighbors, with the purpose of 
avoiding particle pairing (as suggested by Dehnen & Aly 2012). 
The time stepping adopted is described in Durier & Dalla Vecchia 
( 2012 ). 

Gadget 3-X (Murante) 

This project employs two versions of GADGET3-X (G3-X). 
The default code, GADGET3-X-Std (G3X-Std) is a standard ver¬ 
sion of Gadget3 with the cubic spline kernel, 64 neighbours, low- 
order viscosity and a shear flow limiter (Balsara, 1995) and no con¬ 
duction. The improved version GADGET3-X-ArtCond (G3X-Art; 
Beck et al. 2015) is an advanced version of Gadget3 with the 
Wendland kernel functions (Dehnen & Aly, 2012) with 200/295 
neighbours, artificial conductivity to promote fluid mixing follow¬ 
ing Price (2008) and Tricco & Price (2013), but with an additional 
limiter for gravitationally induced pressure gradients, a time-step 
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limiting particle wake-up scheme Pakmor et al. (2012) and a high- 
order scheme for gradient computation (Price, 2012), shock detec¬ 
tion and artificial viscosity similar to Cullen & Dehnen (2010) and 
Hu et al. (2014). A companion paper (Beck et al., 2015) presents 
an improved hydrodynamic al scheme and its performance in a wide 
range of test problems. 

GADGET3-SPHS (Power) 

Gadgets -SPHS (G3-SPHS) was developed to overcome the 
inability of classic SPH to resolve instabilities. It has been imple¬ 
mented in the GADGETS code. The key differences with respect to 
standard GADGETS are in the choice of kernel and in dissipation. 
Rather than the cubic spline kernel, GS-SPHS uses as an alterna¬ 
tive either the HOCT kernel with 442 neighbors or the Wendland 
C4 kernel with 200 neighbors. A higher order dissipation switch 
detects, in advance, when particles are going to converge. If this 
happens, conservative dissipation is switched on for all advected 
fiuid quantities. The dissipation is switched off again once particles 
are no longer converging. This ensures that all fluid quantities are 
single-valued throughout the fiow by construction. 

Gadget3-Magneticum (Saro) 

Magneticum (G3-Magneticum) is based on the entropy- 
conserving formulation of SPH (Springel & Hernquist, 2002). 
Higher order kernel based on the bias-corrected, sixth-order Wend¬ 
land kernel (Dehnen & Aly, 2012) with 295 neighbors. Also in¬ 
cluded is a low viscosity scheme to track turbulence (Dolag et al., 
2005) and isotropic thermal conduction with 1/20^^ Spitzer (Dolag 
et al., 2004). 

Gadget 3-PESPH (February) 

Gad GET 3-PE SPH (G3-PESPH) is an implementation of 
Gadgets with pressure-entropy SPH (Hopkins, 2013) which fea¬ 
tures special galactic wind models. The SPH kernel is an HOCTS 
(n=5) B-spline with 128 neighbors. The time dependent artificial 
viscosity is that of Morris & Monaghan (1997). 

GADGET3-MUSIC (Yepes) 

The original MUSIC runs (G3 -MUSIC) were done with the 
Gadgets code, based on the entropy-conserving formulation of 
SPH (Springel & Hernquist, 2002). GADGETS employs a spline 
kernel (Monaghan & Lattanzio, 1985) and parametrize artificial 
viscosity following the model described by Monaghan (1997). 

GADGET3-OWLS (McCarthy) 

The Gadgets -OWLS (G3-OWLS) simulations were run us¬ 
ing a version of the Lagrangian TreePM-SPH code GADGETS, 
which was significantly modified to include new subgrid physics 
for radiative cooling, star formation, stellar feedback, black hole 
growth and AGN feedback, developed as part of the OWLS/cosmo- 
OWLS projects (Schaye et al., 2010). Standard entropy-conserving 
SPH (Le Brun et al., 2014) was used with 48 neighbors. 

Gadget 2-X (Kay) 

Gadget2-X (G2-X) is a modified version of the Gad- 
GET2 code (SP05), using the TreePM gravity solver and standard 
entropy-conserving SPH. It includes a number of sub-grid mod¬ 
ules to model metal-dependent radiative cooling, star formation and 
feedback from supernovae and AGN. 

Hydra (Thacker) 


Table 2. Key numerical parameters used for each run. Columns 2 and 3 list 
values for the Plummer-equivalent softening lengths for the DM particles 
in physical units; columns 4 and 5 the same but for the gas particles (where 
present); and column 6 the number of EFT cells along each side of the box. 
We use the label ’ Adp’ when adaptive force resolution is used. 


Code name 

cdm 

,max 

^DM 

Cgas 

,max 

‘^gas 

AppT 

CART 

Adp 

Adp 

Adp 

Adp 

256 

Arepo 

33.0 

5.5 

Adp 

Adp 

512 

G2-Anarchy 

20.0 

5.0 

20.0 

5.0 

512 

G3-X 

8.0 

6.0 

8.0 

6.0 

256 

GS-SPHS 

5.0 

5.0 

Adp 

Adp 

1024 

G3-Magneticum 

11.25 

3.75 

3.75 

3.75 

256 

G3-PESPH 

5.0 

5.0 

5.0 

5.0 

256 

G3-MUSIC 

8.0 

6.0 

8.0 

6.0 

512 

G3-OWLS 

9.77 

5.0 

9.77 

5.0 

1024 

G2-X 

24.0 

6.0 

24.0 

6.0 

256 

Hydra 

Adp 

5.0 

Adp 

5.0 

512 


HYDRA-OMP (Hydra) (Thacker & Couchman, 2006) is a 
parallel implementation of the earlier serial HYDRA code (Couch- 
man, Thomas & Pearce, 1995). Aside from the parallel nature of the 
code, HYDRA-OMP differs from the serial implementation by us¬ 
ing the standard pair-wise artificial viscosity along with the Balsara 
modification for rotating flows. Otherwise, the SPH implementa¬ 
tion is ’’classic” in nature, using 52 neighbors, and does not include 
more recently preferred kernels, terms for conduction or explicit 
shock tracking to modify the dissipation. It also uses a conservative 
time-stepping scheme that keeps all particles on the same synchro¬ 
nization. 


2.3 Progresses with modern SPH codes 

Since the first development of SPH by Gingold & Monaghan 
(1977) and Lucy (1977) great advancements have been made to 
improve the accuracy and stability of SPH simulations. In partic¬ 
ular, much attention has been given to the treatment of discon¬ 
tinuities. Artificial viscosity is necessary for a proper fiuid sam¬ 
pling at shocks and to prevent particle interpenetration. The first 
spatially constant low-order formulations of artificial viscosity ap¬ 
plied viscosity not only in shocks, but also in shearing fiows and 
un-shocked regions leading to an over diffusion of kinetic energy. 
Modem formulations of artificial viscosity rely on proper shock de¬ 
tection methods and high-order gradient estimators to distinguish 
between shocked and un-shocked or shearing regions (Morris & 
Monaghan, 1997; Cullen & Dehnen, 2010; Price, 2012; Hu et al., 
2014). Most importantly, they preserve kinetic energy to a much 
higher degree and promote simulations of turbulence or hydrody- 
namical instabilities. Furthermore, SPH intrinsically fails to treat 
different gas phases and their mixing correctly, caused by the lack 
of diffusion terms and an always present spurious surface tension, 
as shown for instance by Agertz et al. (2007). 

Read, Hayfield & Agertz (2010) showed that the mixing prob¬ 
lem in SPH owes to two problems: the force inaccuracy and the lack 
of entropy mixing. Artificial conduction (Price, 2008) or pressure- 
entropy (Ritchie & Thomas, 2001; Saitoh & Makino, 2013; Hop¬ 
kins, 2013) formulations have been developed to overcome these 
issues. They provide for the transport of heat between particles or 
change the basic physical variables. However, in the presence of 
gravitationally induced pressure and temperature gradients, artifi- 
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Table 3. SPH schemata used for the comparison runs. We list the employed kernel functions and number of neighbours (A^sph) as well as the minimum 
smoothing length (/imin) in terms of the gravitational softening length. Furthermore, we provides information about the artificial viscosity and conductivity 
switches. 


Code name 

Hydrodyn. Kernel 

Nsph 

bmin 

Art. Vise. 

Shearflow 

Mixing 

Limiter 

G2-Anarchy 

Wendland C2 

100 ±3 

0 

Adaptive 

LowOrder 

Artificial 

Yes 

G3-XArt 

Wendland C6 

295 ± 10 

0.1 

Adaptive 

HighOrder 

Artificial 

Yes 

G3-SPHS 

Wendland C4 

200 ± 0.5 

1.0 

Adaptive 

LowOrder 

Artificial 

Yes 

G3-Magneticum 

Wendland C6 

295 ± 0.5 

0.001 

Adaptive 

HighOrder 

Physical 

- 

G3-PESPH 

HOCTS B-spline 

128 ±2 

0.1 

Adaptive 

LowOrder 

- 

- 

G3-MUSIC 

Cubic spline 

40 ±3 

0.1 

Constant 

LowOrder 

- 

- 

G3-XStd 

Cubic spline 

64 ± 1 

0.1 

Constant 

LowOrder 

- 

- 

G3-OWLS 

Cubic spline 

48 ± 1 

0.01 

Constant 

LowOrder 

- 

- 

G2-X 

Cubic spline 

50 ±3 

1 

Constant 

LowOrder 

- 

- 

Hydra 

Cubic spline 

53±15 

0.5 

Constant 

LowOrder 

- 

- 


cial conduction schemes might lead to unwanted transport of heat, 
making the use of numerical limiters as well as correction terms are 
highly desirable. Read & Hayfield (2012b) showed that pressure- 
entropy SPH fails for strong shocks and/or if the density gradient is 
large. This was significantly improved by Hopkins (2013) who de¬ 
rived a conservative pressure-entropy SPH for the first time. How¬ 
ever, the problem with large density gradients remained. Read & 
Hayfield (2012b) propose instead to use higher order switching, 
similarly to Cullen & Dehnen (2010), but applied to all advected 
fiuid quantities. This solved the mixing problem in SPH also for 
high density contrast and opened the door to ’’multimass” SPH for 
the first time. 

Lastly, the commonly employed cubic spline kernel function 
can easily become unstable, which leads to the pairing instability, 
incorrect gradient estimators and a poor fiuid sampling. Therefore, 
a change of kernel function is highly recommended, and Wendland 
kernels (Dehnen & Aly, 2012) are now commonly used to retain 
fiuid sampling. Table 3 provides an overview of the different SPH 
codes participating in our cluster comparison project and lists their 
numerical details. 


3 THE SIMULATION 

The cluster we have adopted for this project was drawn from the 
MUSIC-2 sample (Sembolini et al. 2013; Sembolini et al. 2014; 
Biffi et al. 2014) which consists of a mass limited sample ^ of 
re-simulated halos selected from the MultiDark cosmological sim¬ 
ulation (Prada et al., 2012). This simulation is dark-matter only 
and contains 2048^ (almost 9 billions) particles in a (l/i“^Gpc)^ 
cube. It was performed in 2010 using ART (Kravtsov, Klypin & 
Khokhlov, 1997) at the NASA Ames Research centre. All the data 
from this simulation are accessible online through the MultiDark 
Database.^ The run was done using the best-fitting cosmological 
parameters to WMPA7-fBAO-hSNI (Qm = 0.27, Qb = 0.0469, 
Qa = 0.73, as = 0.82, n = 0.95, h = 0.7, Komatsu et al. 2011). 
This is the reference cosmological model used in the rest of the 
paper. 

The MUSIC-2 cluster catalogue was originally constructed by 


^ specifically, it is cluster 19 of MUSIC-2 dataset; all the initial conditions 
of MUSIC clusters are available at http : //music . ft. uam.es 
3 www.MultiDark.org 


selecting all the objects in the simulation box which are more mas¬ 
sive than 10^^ h~^MQ at redshift z = 0. In total, 282 objects 
were found above this mass limit. A zooming technique described 
in Klypin et al. (2001) was used to produce the initial conditions 
for the re-simulations. All particles within a sphere with a radius 
of 6 h~^Mpc around the centre of each selected object at z = 0 
were found in a low-resolution version (256^ particles) of the Mul¬ 
tiDark volume. This set of particles was then mapped back to the 
initial conditions to identify the Lagrangian region corresponding 
to a 6 /i“^Mpc radius sphere centered on the cluster centre of mass 
at z = 0. The initial conditions of the original simulations were 
generated on a finer mesh of size 4096^. By doing so, the mass res¬ 
olution of the re-simulated objects was improved by a factor of 8 
with respect to the original simulations. In the high resolution re¬ 
gion the mass resolution for the dark matter only simulations corre¬ 
sponds to mDM=1.09x 10^ /i^^Mq. For the runs with gas physics, 
mDM=9.01xl0® h~^MQ and to mgas=L9x 10® /i^^M©. 


4 DARK MATTER RUN COMPARISON 

We first completed a dark matter only simulation, performed using 
the parameter values given in Table 2. These were chosen indepen¬ 
dently by each modeling group following their previous experience. 
A visual comparison of the density field centered on the cluster at 
z = 0 is given in Figure 1 . While it is clear that all the methods 
have followed the formation of the same object (with a significant 
improvement with respect to Figure 1 of SB99) the precise loca¬ 
tion of the major subhalo is not accurately recovered. For several 
methods it has already crossed R 2 oo^ (the radius enclosing an over¬ 
density of 200 relative to the critical density) while in others it is 
still falling in. The variance across this figure illustrates the typical 
range of outcomes from commonly used cosmological codes. The 
major cause of the discrepancy (for the Gadget based codes at 
least) is the size of the base level particle mesh. Those implemen¬ 
tations which employed a base mesh of 256® did not sufficiently 
resolve the interface region between this low resolution mesh and 
the higher resolution region placed over the cluster. Improving the 
resolution of the base level mesh to 512® alleviates this problem 
and aligns the dark matter component well. We demonstrate this by 
showing the effect of changing the size of the particle mesh for the 
G3-MUSIC code in the appendix (Figure 13) together with a simi¬ 
lar set of visual images taken from the following non-radiative sim¬ 
ulation when the simulation parameters infiuencing the accuracy of 
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the dark matter are aligned across the code. This demonstrates that, 
given a common set of parameters, the dark matter framework un¬ 
derlying the simulation can be made to look very similar (this is 
not surprising at least for the Gadget based codes and Arepo as 
they all use the same gravity solver). 

Figure 2 displays the radially averaged projected dark mat¬ 
ter density profiles for the 12 different non-aligned DM-only runs, 
including also the residuals relative to the density profile of the 
reference G3-MUSIC simulation. The secondary peak marks the 
location of the major subhalo, at ~ l/i“^Mpc, significantly 
closer to the centre in some simulations due to the size of the top 
level particle mesh employed. All the simulations except CART 
lie well within 10 per cent of each other at all radii with the Hy¬ 
dra simulation being indistinguishable from the Gadget runs. 
CART produces a cluster that is slightly more centrally concen¬ 
trated than for the particle based approaches, especially within the 
inner 100/i“^kpc. 

The subhalo mass function at z = 0 (Figure 3) is 

recovered with a very close agreement (differences are al¬ 
ways below 20 per cent at all masses) by all codes. Subha¬ 
los here have been identified using AHF (Gill, Knebe & Gib¬ 
son 2004; Knollmann & Knebe 2009; freely available from 
http://popia.ft.uam.es/AHF). The number of subhalos 
is essentially identical except for tiny mass differences which are 
driven by the small divergences in radial location that were iden¬ 
tified above. These code to code variations lead to differences in 
the mass associated with each subhalo as the threshold that defines 
where a subhalo is separated from the background halo varies. As 
expected, subhalos closer to the centre of the main halo than the 
equivalent subhalo in one of the other models have a lower recov¬ 
ered mass. 

Comparison of the dark matter distribution and of its radial 
density profile at z = 1 give results similar to those described 
above. We conclude that the typical range of chosen parameters 
for cosmological simulation codes introduces a variation of around 
10 per cent in the density profile of collapsed objects. This scatter 
can be reduced to less than 5 per cent by aligning the numerical ac¬ 
curacy of our calculations. Although this is not essential for many 
applications, we choose to do this in the remainder of this paper so 
that, as we show in the appendix, the underlying dark matter frame¬ 
work agrees closely, allowing us to focus on differences resulting 
from the various hydrodynamical schemes. 


5 NON-RADIATIVE RUN COMPARISON 


We now proceed to include a gas phase into our calculations. We 
repeat the simulation of the same cluster as used in Section 4 in¬ 
cluding gas which however we do not allow to radiate energy. Fig¬ 
ure 4 shows some of the global properties of the selected cluster as 
calculated by the different codes used in this work: radius, mass, 
mass-weighted temperature, gas fraction, dark matter velocity dis¬ 
persion and axial ratios. All these quantities have been calculated 
at an aperture radius corresponding to , the radius enclosing 
an overdensity of 200 relative to the critical density, defined as 


Pc{z) 


SHjEjzf 

SttG 


( 2 ) 


where Hq is the present value of the Hubble constant, G is the uni¬ 
versal gravitational constant - using the same definition we refer in 
the text to i^ 25 oo and Rioo^ as the radii enclosing an overdensity 
of 2500 and 500 to Pc(z). It is interesting to note that all the codes 


were able to recover the same values for the different properties of 
the halo with a scatter smaller than 1 per cent for mass, radius, axial 
ratio and dark matter velocity dispersion and around 2 per cent for 
baryon fraction and gas temperature. The same properties were es¬ 
timated with a scatter between 5 and 10 per cent in SB99, with dif¬ 
ferences of up to 30 per cent between the maximum and minimum 
values: in our comparison the same difference is always below 5 
per cent (only for the gas fraction we register a disagreement of 8 
per cent between the maximum and minimum values). 

Thumbnail images of the gas density distribution for each of 
the methods at z = 0 are given in Figure 5. We see a dramatic 
variation in the central concentration of the gas, with some methods 
having significantly larger extended nuclear gas regions. 

This trend is born out by the radial gas density profiles given 
in Figure 6. We see that the radial gas density is significantly more 
extended for CART and Arepo when compared to the traditional 
SPH schemes employed by some SPH codes such as Hydra and 
G2-X. SPH codes that implement artificial diffusion are quite close 
to CART and AREPO. Lagrangian methods with entropy mixing 
(including AREPO which is Lagrangian in spirit) are always very 
close to each other, while CART produces a gas density profile that 
is shallower in the central regions and steeper at larger radii. 

The difference in the radial gas density compared to the fidu¬ 
cial G3-MUSIC run are shown in the top panel of Eigure 6. All 
the residuals are calculated to the density profile of the reference 
G3-MUSIC simulation (and we adopt this definition for all the ra¬ 
diative profiles shown from now on). At z = 0 the lowest central 
densities are an order of magnitude lower than in the G3-MUSIC 
simulation while the highest central densities are around 5 times 
larger, i.e. the variation in the central gas density across our runs is 
nearly two orders of magnitude. The scatter becomes more moder¬ 
ate when considering the outer region of the cluster, not exceeding 
20 per cent at radii larger than R^boo- 

We next show the radial temperature profiles for all the simu¬ 
lations in Eigure 7. We use the mass-weighted temperature, defined 
as: 


_ '^iTirrii 

T,i^i 


(3) 


where rrii and Ti are the mass and the electronic temperature of 
the gas particles . The central temperatures vary by more than a 
factor of 3, with a group of methods displaying a central tempera¬ 
ture inversion with inner temperatures around half the peak value 
of 7-8 keV. In contrast, some codes display a monotonically rising 
temperature profile as the radius falls with a peak temperature up 
to 14 keV at the very centre. This group of codes consists of those 
with the most extended radial gas density profiles. Also in this case 
CART results are different from AREPO and SPH with thermal 
diffusion. It is interesting that the jump in temperature seen at ~1 
h~^Mpc for all Lagrangian methods is not seen by CART, proba¬ 
bly because in this case the substructure that causes this bump is 
located in a slightly different halo position due to the differences in 
the merger phase. At radii larger than R^lm the scatter is signifi¬ 
cantly more moderate, and the residuals appear to be a factor of 2 
smaller than in Eigure 17 of SB99 in the same cluster regions. 

Einally we show the radial gas entropy profiles for all the 
codes in Eigure 8. We adopt the definition of entropy commonly 
used in the literature by works on X-ray observations: 


S{R) 


Tgas{R) 

nV\R) 


(4) 


where Ue is the number density of free electrons of the gas. This 
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Figure 1. Projected density of the dark matter halo at 2 ; = 0 for each simulation as indicated. The box is 2h ^Mpc on a side. The white circle indicates 
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Figure 2. Radial density profiles for the dark matter only simulations at 2 ; = 0 (bottom panel) and difference between the radial density profiles of each dark 
matter only simulations at 2 ; = 0 and the reference G3-MUSIC density profile (top panel). The vertical dashed line corresponds to R 2500 and vertical dotted 
line to i? 5 oo of the reference G3-MUSIC values. 


displays the now classic split between grid-based codes and tra¬ 
ditional SPH methods, such as Hydra and G2-X, which show a 
falling inner entropy as the radius is decreased all the way into the 
very centre. This is completely consistent with the inner tempera¬ 
ture inversion and high central density. Conversely, the grid-based 
codes such as CART and Arepo display the well known flat in¬ 
ner entropy cores that result from rising inner temperature profiles 
and extended gas densities. However, we see that in between these 
extremes we have a full range of entropy profiles that depend on 
the specific SPH implementation employed. Differently from what 
was shown in SB99 (see their Figure 18) modern, sophisticated 
SPH codes which employ mixing are now capable of recovering 
entropy profiles that lie anywhere between the core-less traditional 
SPH schemes and the cored profiles of the grid-based approaches 
depending upon the precise nature of the scheme and the amount 
of mixing employed. We highlight that modern SPH codes such as 
G3-SPHS, G2-Anarchy and G3-XArt are able to recover the same 
flat entropy core observed for CART and Arepo with a scatter 
smaller than 20 per cent, even in the inner cluster regions. G3- 
PESPH and G3-Magneticum, which have artificial viscosity switch 
but different artificial conductivity with respect to the other modern 
SPH codes, show and intermediate behavior between standard and 
modern SPH codes. 


5.1 Other quantities in the non-radiative simulations 


It is important to note that the differences in radial gas density, 
temperature and entropy evidenced above are not driven by code 
issues such as poor thermalization or large scale flows. In Figure 9 
we show the ratio of gas thermal U to kinetic energy K Sit z = 0: 


radial profile of all the simulations. All the methods agree closely 
on the value of 77 as a function of halo radius and none display any 
evidence of poor thermalization. Interestingly, we note that CART 
is the most efficient in converting kinetic energy into thermal en¬ 
ergy (though also the different merger phase may contribute to this 
offset). Even considering the difference between CART and the 
other codes, the scatter is always below 20 per cent. 

Given our radial dark matter density and gas density profiles 
we can also calculate the radial gas fraction for all the methods. 
In Eigure 10 we show the radial profiles of the depletion factor T, 
defined as: 


Mga,{< R) 

M{<R) \Q.m) 


(6) 


The results reflects the existence of methods that produce very cen¬ 
trally concentrated gas and methods yielding an extended core with 
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Figure 3. Subhalo mass function for the dark matter only simulations at 2; = 0 (bottom panel) and difference between the subhalo mass function of each dark 
matter only simulations at 2; = 0 and the reference G3-MUSIC subhalo count (top panel). 


much higher average entropy. For those codes with an extended en¬ 
tropy core the gas fraction falls significantly with decreasing radius 
and can reach values as low as 5 per cent of the Universal baryon 
fraction for these non-radiative simulations. Within 100/i“^kpc it 
can fall below 25 per cent of the cosmic gas fraction in both grid- 
based and modem SPH codes. The differences in the radial gas 
fraction reflect those detected in the gas density profiles and warn 
about using a universal calibration of the baryon depletion at i? 25 oo 
based on simulations in cosmological applications of the cluster 
baryon fraction, especially when using only non-radiative gas. We 
expect these results to be very different (and closer to observational 
results) when radiative physics is included. The scatter in the outer 
regions of the cluster appear by the way to be much more limited 
(less than 20 per cent) with respect to the results shown in Figure 
13 of SB99, where differences of up to 50 per cent between the 
different codes where registered even close to the virial radius. 

We can combine our measurements of the gas density and 
temperature to produce Figure 11, the gas pressure profiles and 
Figure 12, the X-ray emission profiles. We define the pressure as 
P = PgasT and we normalize the profiles to the value of P500 
(the value of the pressure as calculated at R^qq) in order to be con¬ 
sistent with the definition of universal pressure profile introduced 
by Amaud et al. (2010). The X-ray luminosity profile is defined 
as 4:7tR^ Lx and we approximate the X-ray luminosity density as 
Lx = The variation in the gas density and temperature 


produce very different pressure and X-ray emission profiles. As ex¬ 
pected, the pressure profiles continue to rise all the way into the 
centre for all codes (as the central gas is close to hydrostatic equi¬ 
librium in all cases). Due to the very high central density, the central 
X-ray emission for Hydra is over two orders of magnitude higher 
than that found by the grid-based and some modern SPH meth¬ 
ods which form extended cores. The differences observed in CART 
profiles, especially in Figure 11, can be attributed to the different 
merger phase. 

6 CONVERGENCE AND SCATTER BETWEEN 
SIMULATIONS 

6.1 Dark matter only runs 

The major cause of the discrepancy (for the Gadget based codes 
at least) is the size of the base level particle mesh. A base mesh 
of 256^ is not sufficient to resolve the interface region between the 
low resolution region (i.e. the base level) and the high resolution 
region (i.e. the refined level) placed over the cluster; we found that 
a base level of at least 512^ is required to ensure that the dark mat¬ 
ter component is well aligned across the different codes. For the 
N -body only simulations we find excellent agreement between the 
density profiles of the main halo and the statistics of their subha¬ 
los (Figures 2 and 3), once the input parameters of the various runs 
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Figure 4. Global properties of the cluster calculated by the different codes. All quantities are computed within i^ 200 ^- From top panel to bottom panel: (1) the 
one-dimensional dispersion of the dark matter, (2) the axial ratio (h/a in black, cjam red), (3) the mass-weighted temperature, (4) the gas fraction (dotted 
line represented the value of the cosmic ratio from WMAP7 (Komatsu et ah, 201 1), (5) the radius and (6) the total cluster mass. The solid lines represent the 
median value for each one of the plotted quantities and the dashed lines the -i-/- 1 per cent scatter. 
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Figure 5. Projected density of the gas halo at 2 ; = 0 for each simulation as indicated. The box is 2/i“^Mpc on a side. The white circle indicates for 

the halo, the black circle shows the same but for the G3-MUSIC simulation. 
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R [h ' Mpc] 


Figure 6. Radial gas density profiles at 2 ; = 0 for each simulation as indicated (bottom panel) and difference in radial gas density profiles at 2 ; = 0 between 
each simulation and the reference G3-MUSIC simulation . The vertical dashed line corresponds to i ^2500 and the dotted line to i^soo of the reference G3- 
MUSIC values. The error bars on G3-SPHS (black) and G3-MUSIC (red) are calculated from the scatter between snapshots averaged over 0.27 Gyrs. The data 
are cut off when the radius goes below the softening scale of the code. 


are matched (see Appendix). However, even the ’matched’ simula¬ 
tions show important phase differences between simulations (Fig¬ 
ure Al). The numerous implementations of Gadget agree remark¬ 
ably well, but this is perhaps not surprising given that the gravity 
solver is in each case built on the same foundation. CART, how¬ 
ever, shows a marked disagreement, as a large sub halo appears to 
be absent. Similar results were reported recently in Power, Read & 
Hobbs (2014) when comparing G3-SPHS with the Ramses AMR 
code (see their Figure 7). It is not clear that this is either a problem 
or unexpected. We are, after all, modeling a chaotic system and 
while we can expect the statistics of the density field to agree be¬ 
tween codes, it is likely unreasonable to expect the precise merger 
history and merger phase to agree between codes for a single ob¬ 
ject. Ideally, we would compare the halo statistics between codes 
for a large ensemble of cluster zoom simulations. Such a study is 


beyond the scope of this present work and will treated in an future 
paper dedicated to the study of the subhalos. 


6.2 Adding gas: AMR versus SPH 

Issues such as the merger phase appear to affect the resulting dark 
matter properties at the ~ 10 per cent level (Figures 2-3). However, 
once adding gas, the merger phase produces a larger effect than this 
on the resulting gas distribution. This is because entropy in clusters 
is generated in shocks during collapse and mergers (Mitchell et al., 
2009; Power, Read & Hobbs, 2014). The different codes explored 
here broadly fall into four categories. There are the ‘classic’ SPH 
methods (Hydra; G3-MUSIC; G3-OWLS; G3-XStd; G2-X); the 
’modern’ SPH methods that attempt to correct for problems with 
mixing in classic SPH (G3-XArt; G3-SPHS; G2-Anarchy; G3- 
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Figure 7. Radial temperature profile at 2 ; = 0 for each simulation as indicated and difference between each simulation and the reference G3-MUSIC simulation. 
The vertical dashed line corresponds to i ^2500 and the vertical dotted line to R 500 of the reference G3-MUSIC values. 


PESPH; G3-Magneticum); an Adaptive Mesh refinement method 
(CART) and a hybrid ‘moving mesh’ method (Arepo; see section 
2 for further details of each method and the differences between 
them). The AMR code CART that is the most discrepant in its dark 
matter distribution is also the most discrepant in in its gas distribu¬ 
tion, particularly at large radii. This is almost certainly attributable 
to the different merger phase, as can be seen from Figure Al. 

6.3 The trouble with ‘classic’ SPH 

Despite the excellent agreement between the ‘dark matter only’ 
runs, the classic SPH simulations also show scatter in their central 
gas density, temperature and entropy that is significantly larger than 
the theoretical error on each (as calculated from the scatter between 
snapshots averaged over 0.27 Gyrs; see the error bars marked on the 
plots shown in Figures 6 and 8). This is likely because in classic 
SPH, low entropy particles artificially sink to the cluster centre due 
to the lack of mixing (Power, Read & Hobbs, 2014), amplifying the 
effect of small changes in merger phase. We see this in Figure 5. 

6.4 Modern SPH and Arepo 

Of considerable interest are the differences between the modern 
SPH flavors (G3-XArt; G3-SPHS; G2-Anarchy; G3-PESPH; G3- 
Magneticum) and the moving mesh code Arepo. There is excel¬ 


lent agreement between the dark matter distributions across the 
runs, which allows us to isolate the effect of the hydrodynamics 
solver (Figure Al). For the gas distributions shown in Figure 6, 
Arepo G3-XArt and G3-SPHS are in excellent agreement with 
one another, agreeing almost perfectly within our estimated the¬ 
oretical error, while G2-Anarchy seems to be the outlier. In Fig¬ 
ure 7 the temperature profiles for Arepo and G3-XArt are very 
closely matched but look to have an offset from G3-SPHS and G2- 
Anarchy, which are very close to one other. In Figure 8 the radial 
entropy profiles of G3-SPHS, G3-XArt and Arepo are again very 
close to one another while G2-Anarchy is more discrepant, though 
closer to CART. The origin of these differences is yet to be ex¬ 
plained, although we note that it cannot be attributed to different 
merger phases and must result from the hydrodynamics solver. In 
the case of G2-Anarchy a possible cause of discrepancy may be the 
choice of the kernel (using a C4 kernel with 200 neighbors gives 
slightly different values for the central entropy). Power, Read & 
Hobbs (2014) showed that at the resolution of the simulations in 
this paper, G3-SPHS is numerically converged. It would be inter¬ 
esting to see if the differences between G3-SPHS/G3-XArt/Arepo 
and G2-Anarchy remain with increasing resolution. We defer such 
a resolution study to future work, whose intent will be to narrow 
down these numerical uncertainties. 

In the Arepo simulations we have used the total energy as 
a conserved quantity in the Godunov scheme, which is the usual 
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Figure 8. Radial entropy at 2 ; = 0 (bottom panel) for each simulation as indicated and difference between each simulation and the reference G3-MUSIC 
simulation (top panel). The dashed line corresponds to i^2500 and the dotted line to -R500 of the reference G3-MUSIC values. The error bars on G3-SPHS 
(black) and G3-MUSIC (red) are calculated from the scatter between snapshots averaged over 0.27 Gyrs. 
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Figure 9. The ratio of kinetic to thermal energy in the gas, 77 , measured radially at 2 ; = 0 for each simulation as indicated (top) and difference between each 
simulation and the reference G3-MUSIC simulation (bottom). Dashed line corresponds to R 2500 and dotted line to R 500 of the reference G3-MUSIC values. 


choice in finite volume codes and has been used in most recent 
Arepo studies (as, e.g., in Vogelsberger et al. 2014). As discussed 
in Springel (2010), using an entropy-energy formalism results in 
smaller entropy cores and higher central gas densities somewhat 
closer to classic SPH results (similar results are also shown by 
GIZMO, Hopkins 2014). Following most recent Arepo studies, we 
have not employed the later method here due to concerns regarding 
the artificial suppression of weak shocks and the potentially less 
accurate energy conservation. 

Finally, we note that the results of G3-PESPH are very dif¬ 
ferent from those of the other modern SPH flavours (with the 
exception of G3-Magneticum), and are more similar to those of 
classic SPH simulations. A key difference is that this version of 
G3-PESPH does not include any explicit conductivity or mixing, 
while all the other modern variants do. Hu et al. (2014) showed 
that PESPH performs much better than previous versions of SPH 
for surface instabilities by greatly mitigating surface tension prob¬ 
lems, but in areas of very strong shocks (M~1000) adding artificial 
conduction provides a better match to analytic solutions. Insights 
into the behavior of G3-PESPH may be gained by considering the 
OSPH method presented in Read, Hay field & Agertz (2010), and 
the earlier multiphase RT method by Ritchie & Thomas (2001). As 
pointed out by Read & Hayfield (2012a), the RT method only works 
well for relatively small entropy contrasts between different fluid 
phases. As the entropy contrast becomes very large, the admixture 
of low and high entropy particles within the smoothing kernel cre¬ 
ates large pressure fluctuations that prevent mixing, as in classic 
SPH. This was recognized also by Ritchie & Thomas (2001) who 
proposed scaling the neighbor number with the entropy contrast to 
combat this. However, this rapidly becomes prohibitively expen¬ 


sive in realworld applications, which led Read & Hayfield (2012a) 
to abandon such RT methods in favor of dissipation switching, as 
proposed by Price (2008); such switching is common to G3-XArt, 
G3-SPHS and G2-Anarchy and helps to explain their similarity. 
The discrepancy between G3-PESPH and the other modern SPH 
flavours may also reflect the treatment of artificial viscosity adopt¬ 
ing the artificial viscosity model suggested by Cullen & Dehnen 
(2010) can produce larger entropy gains in shocks, but the authors 
of G3-PESPH do not use this because it also seems to add substan¬ 
tial entropy into very diffuse intergalactic gas that may be spurious 
(N. Katz, priv. comm.). In short, it is unclear how much artificial 
conductivity and/or mixing is appropriate in SPH, or even whether 
the mesh codes are providing the correct solution that all SPH codes 
should be trying to emulate. Nonetheless, consistency with mesh 
codes appears to require modern SPH using conduction, mixing, 
and/or a higher-order dissipation switch. 

The discrepancy between G3-PESPH and the other modem 
SPH flavours may also reflect the treatment of artificial viscosity 
- adopting the artificial viscosity model suggested by Cullen & 
Dehnen (2010) should lead to a better agreement with mesh-based 
codes. 


7 SUMMARY & CONCLUSIONS 

We have investigated the performance of 12 modem astrophysi- 
cal simulation codes - CART, Hydra, Arepo and 9 versions of 
Gadget with different SPH implementations - by carrying out 
cosmological zoom simulations of a single massive galaxy cluster. 
Our goal was to assess the consistency of the different codes in re- 
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Figure 10. Radial gas fraction at 2 ; = 0 relative to the cosmic value for each simulation as indicated (top) and difference between each simulation and the 
reference G3-MUSIC simulation (bottom). The dashed vertical line corresponds to R 2500 and dotted vertical line to R 500 of the reference G3-MUSIC values. 


producing the spatial and thermodynamical structure of dark matter 
and non-radiative gas in the cluster. 

As our initial step, we ran dark matter only versions of the 
simulations with each code using its preferred set of numerical pa¬ 
rameters (e.g. time step accuracy, gravitational softening, dimen¬ 
sion of the particle mesh), and examined the spherically averaged 
mass density profile and the spatial distribution of substructures. 
We found good consistency between the density profiles recovered 
by the codes at approximately the 10 per cent level, while there 
were small variations in the positions of substructures. When these 
simulations were re-run with a common set of numerical parame¬ 
ters, we found that these small variations could be suppressed (es¬ 
sentially entirely, in the case of the Gadget codes). 

By adopting this common parameter set, we were able to 
isolate those differences between the results of the hydrodynami- 
cal simulations that arise only from the choice of hydrodynamical 
solver, rather than from the complex interplay of the hydrodynam¬ 
ical and gravity solvers. Interestingly, we found that the resulting 
gas density profiles varied substantially amongst the codes. Our key 
findings can be summarized as follows: 

• Some codes, essentially the oldest, with classic SPH imple¬ 
mentations, exhibit continually falling inner entropy profiles, with¬ 
out any evidence of an entropy core. This is because these codes, 
particularly Hydra ,were carefully designed to be entropy con¬ 
serving with very low levels of mixing. This lack of mixing pre¬ 
serves low-entropy gas particles at the centers of all objects, in¬ 
cluding subhalos, which survive until late times. As the cluster re¬ 
laxes, these particles sink to the centre of the radial density profile, 
decreasing the central entropy. 


• In contrast, the grid-based codes CART and Arepo produce 
extended cores with a large constant entropy core. In these mesh 
based codes mixing of entropy arises as a consequence of the nu¬ 
merical diffusion associated with the Riemann solver: they natu¬ 
rally mix entropy between gas elements, essentially eliminating the 
very low entropy material. 

• Modern SPH codes such as G3-ANARCHY, G3-SPHS and 
G3-XART, which have dissipative switches and new kernels, can 
bridge the gap between the classic SPH codes and grid based codes, 
and produce entropy cores that are indistinguishable from those of 
the grid-based codes. 

Our results confirm that the discrepancies between grid-based 
codes and SPH codes in describing the radial entropy profile of sim¬ 
ulated clusters, identified by the Santa Barbara comparison project 
presented in Frenk et al. (1999), can be overcome by modern SPH 
codes. Importantly, all the codes employed in this work succeed in 
recovering the global properties and most of the radial profiles of a 
simulated large galaxy cluster with much greater accuracy and sig¬ 
nificantly smaller scatter than those presented in Frenk et al. (1999); 
this highlights the enormous strides in the development of astro- 
physical hydrodynamical simulation codes over the last decade. 

This work constitutes the first in a series of papers in which we 
examine in detail the predictions of modern astrophysical hydrody¬ 
namical simulation codes. The next paper in this series will focus 
on simulations of the same galaxy cluster, now modeled with a va¬ 
riety of galaxy formation processes including cooling, star forma¬ 
tion, supernovae, and feedback from active galactic nuclei (AGN). 
This will allow us to establish how radiative processes affect the 
entropy cores of simulated clusters. Subsequent papers will look at 
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Figure 11. Radial gas pressure at 2 ; = 0 measured in each simulation as indicated (bottom panel) and difference between each simulation and the reference 
G3-MUSIC simulation (top panel). The pressure, as well as the radius, has been normalized to the corresponding value at -R 500 for each code. The dashed 
vertical line corresponds to R 2500 of the reference G3-MUSIC value. 


the recovery of cluster properties such as X-ray temperature and 
Sunyaev-ZeTdovich profiles; gravitational lensing; and cluster out¬ 
skirts and hydrostatic mass bias; all of which will add to our un¬ 
derstanding of how consistently the results of different codes can 
inform our understanding of galaxy cluster properties. 
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Figure 12. Radial X-ray luminosity profiles at 2 ; = 0 for each simulation as indicated (bottom panel) and difference between each simulation and the reference 
G3-MUSIC simulation (top panel)..The dashed vertical line corresponds to R 2500 of the reference G3-MUSIC value. 
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8 DARK MATTER ALIGNMENT 

In order to perform a clean comparison of the various gas physics 
implementations we carefully aligned the underlying gravitational 
framework for each model. While Eigure 1 illustrates the range of 
outcomes that result from a blind comparison using individual pa¬ 
rameter choices, we can choose a common parameter set for those 
quantities that control the accuracy of the gravitational forces. Eor 
instance, for Gadget, Table 4 gives the parameter choices made 
independently by each group. The final row lists the common pa¬ 
rameter set adopted for the non-radiative comparison. Eor this com¬ 
mon choice the gravitational evolution of the nine Gadget simu¬ 
lations and Arepo is, as expected, essentially indistinguishable, as 
illustrated by Eigure 13. Eigure 14 shows the radial density dis¬ 
tribution and the difference relative to the G3-MUSIC simulation. 
Eor Hydra, the central gas density is so high that it steepens the 
central dark matter distribution relative to the other codes. 
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Table 4. Numerical parameters used for the Arepo and Gadget runs: accuracy of the time step (ETIA), time step displacement factor (MRDF), minimum 
(MINT) and maximum (MAXT) time step, tolerance on the force accuracy (ETFA), accuracy of the tree algorithm (ETT), Courant factor (CFAC) and double 
precision (DP, DF). 


Code name 

ETIA 

MRDF 

MINT 

MAXT 

ETFA 

ETT 

CFAC 

DP 

DF 

Arepo 

0.025 

0.125 

0 

0.01 

0.0025 

0.6 

0.3 

Y 

Y 

Gadget2-X 

0.02 

0.25 

10-'^ 

0.025 

0.0025 

0.3 

0.15 

Y 

Y 

Gadget3-Magneticum 

0.05 

0.25 

0 

0.05 

0.005 

0.45 

0.15 

Y 

Y 

Gadget3-MUSIC 

0.01 

0.5 

0 

0.01 

0.01 

0.4 

0.15 

Y 

Y 

Gadget3-OWLS 

0.025 

0.25 

0 

1 

0 

0.025 

0.005 

0.6 

0.15 

Y 

N 

Gadget3-SPHS 

0.03 

0.5 

0 

0.02 

0.005 

0.5 

0.4 

N 

N 

Gadget3-X 

0.01 

0.5 
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0.01 

0.01 

0.45 

0.15 

N 
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Gadget3-PESPH 
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Y 

Gadget2-Anarchy 

0.01 
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0.0 
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0.025 
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Common parameter set 
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Figure 14. Radial density prohles for the non radiative simulations at 2 ; = 0 (bottom panel) and difference between the radial density profiles of each non 
radiative simulations at 2 ; = 0 and the reference G3-MUSIC density profile (top panel). The vertical dashed line corresponds to i ^2500 and the vertical dotted 
line to i? 5 oo of the reference G3-MUSIC values. 































